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Quantitative analysis of clumps in the tidal tails of 
star clusters 
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ABSTRACT 

Tidal tails of star clusters are not homogeneous but show well defined clumps in 
observations as well as in numerical simulations. Recently an epicyclic theory for the 
formation of these clumps was presented. A quantitative analysis was still missing. We 
present a quantitative derivation of the angular momentum and energy distribution 
of escaping stars from a star cluster in the tidal field of the Milky Way and derive the 
connection to the position and width of the clumps. For the numerical realization we 
use star- by-star ^V-body simulations. We find a very good agreement of theory and 
models. We show that the radial offset of the tidal arms scales with the tidal radius, 
which is a function of cluster mass and the rotation curve at the cluster orbit. The 
mean radial offset is 2.77 times the tidal radius in the outer disc. Near the Galactic 
centre the circumstances are more complicated, but to lowest order the theory still 
applies. We have also measured the Jacobi energy distribution of bound stars and 
showed that there is a large fraction of stars (about 35%) above the critical Jacobi 
energy at all times, which can potentially leave the cluster. This is a hint that the 
mass loss is dominated by a self-regulating process of increasing Jacobi energy due to 
the weakening of the potential well of the star cluster, which is induced by the mass 
loss itself. 
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1 INTRODUCTION 

Recently well-defined clumps were obse r ved in the tidal 
tails of globular clusters by Leon et al. {2 000) for NGC 



6254 and Pal 12 and by lOdenkirchen et al. I (|200ll . 120031 ) 



for Pal 5. External perturbations like crossing of the galac- 
tic disc, peri-centre passage or near-encounters with other 
globular clusters were discusse d as the source of the clumps 
Capuzzo Dolcetta et al\ 120051 ) . In ICapuzzo Dolcetta et al\ 



2005;) the formation of these clumps in tidal tails of star 



clusters on eccentric orbits were confirmed by numerical 
studies. But the clumps occur also in the tidal tails of star 
clusters on circular orbits with no external push. Recently, 
iKiipper et al\ (|2008l ) presented a theoretical explanation for 
the clump formation in a constant tidal field. It is essen- 
tially due to the epicyclic motion of the stars lost by the 
star cluster. 

We present a quantitative analysis of the tidal tail struc- 
ture for star clusters moving on a circular orbit in the galac- 
tic disc. The analysis is based on numerical simulations with 
realistic particle numbers including an initial mass function 



(IMF) and stellar evolution. We compare the results for star 
clusters at the solar circle and near the Galactic centre. Since 
the epicycle theory is a perturbation theory with respect to 
a circular orbit with constant tidal field, we cannot apply it 
for predictions of clump distances to the eccentric orbits of 
the observed globular clusters. On the other hand the ob- 
servation of tidal tail clumps of open clusters in the galactic 
disc are hampered by the overwhelming number of field stars 
with similar properties. For an identification the contrast 
in density and velocity with respect to the field stars may 
be too small. Additionally the tidal tails may be destroyed 
quickly by the same gravitational scattering process, which 
is also responsible for the dynamical heating of the stellar 
disc. We discuss the observability further in Sect. [5] 

In Section[2]we present the epicyclic theory for the stars 
in the tidal tails and the connection to the mass loss and 
the orbit of the star cluster. In Section [3] we present the nu- 
merical codes used and the properties of the star clusters. 
Section [4] contains the quantitative comparison of the nu- 
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Figure 1. Visualization of the local cartesian (top) and local 
cylindrical coordinates (bottom) for the snapshot at t=1.414Gyr 
of model 10, where the origin is moved from the galactic centre 
to the star cluster centre. Black (blue) dots are the bound stars, 
dark grey dots to the left (red) are the stars of the leading arm 
and light grey dots to the right (green) are the stars in the trailing 
arm. The light dashed (cyan) circle marks the tidal radius. 



merical results with the theoretical predictions. In Section [5] 
we summarize our results. 



2 DYNAMICS OF ESCAPING STARS 

In the first part of this section we derive the orbital prop- 
erties of the stars in the tidal tails in terms of angular mo- 
mentum and energy. Then we discuss the connection of mass 
loss with the Jacobi energy distribution of the stars. Finally 
the tidal tail properties are determined with respect to the 
star cluster orbit. 

Classically the gravitational potential <fr g (R) and ki- 
netic energy Q 2 R 2 /2 are approximated to second order in 
the variable r = R — Ro, which we will also use. For the 
determination of the tidal tail properties with respect to 
the star cluster orbit in Section 12.31 we switch to the Tay- 
lor expansion with respect to angular momentum (e.g. for 
Ro{L) — \/L/Q,(L) equation IA9[) . because L is a constant 
of motion for the tidal tail stars and easily measurable. 



2.1 Motion of tidal tail stars 

As soon as the gravitational potential of the star cluster is 
negligible, stars in the tidal tails move in the axisymmetric 
potential of the Galaxy. The orbits can be calcula ted in the 
epicyclic approximation (see iKiipper et al\ (|2008h for a re- 
cent application to tidal tail structure). The radial offset of 
the epicyclic motion relative to the orbit of the star cluster is 
first order in the angular momentum. The radial amplitude 
depends on the energy excess, which is of second order. 

Since the tidal tails may extend over a considerable 
range in azimuth, we use polar coordinates R, tp, z with ori- 
gin at the galactic centre. Fig.[T]shows the effect of switching 
from local cartesian (awt, 2/cart) to local polar (x,y) coordi- 
nates. We restrict the investigation to orbits in the galactic 
plane z — 0, but a generalization is straightforward. Energy 



E and the z-component of angular momentum L = L z are 
isolating integrals of motion in the axisymmetric potential 
of the Galaxy. The motion is a 2-dimensional harmonic os- 
cillation with the epicyclic frequency k (see equation IA1I) . 
We use the normalized epicyclic frequency f3 = k/Q. The 
'epicentre' (guiding centre) of the oscillation is described by 



(R,<p)a = (Ro,n t) 



(1) 



where Q(R) is the angular frequency of the galactic rotation 
and f2o = fl(Ro)- Ro is determined by the angular momen- 
tum of the star via 



L — Lo = f^o^o 



(2) 



In one epicyclic period T = 2% / no the epicentre moves along 
the circle with radius Ro by 



2-7T 2-7T 

D (T) = — QoRo = -rRo 

K Po 



(3) 



The energy of the star determines the amplitude of the 
oscillation. At the apo- and pericentre Rm, where the radial 
motion vanishes, it can be written as 



E = $ g (R n 



2RI 



(4) 



with the galactic potential <& S (R). Relative to the circular 
motion with 



E = $ s (Ro) + ^ 



(5) 



the radial amplitude r„ 
energy excess 



Rm — Ro is determined by the 



f2 



AE = E-E = %(R m ) - $ g (_R ) + ^ ( — - ) (6) 



ir- 



To second order in r (see App. [A} we find for the radial 
amplitude r m and peri/apocentre position R m 



Rn 



Ro ± r m = Ro ± 



V2AE 
Potto 



(7) 



Note that the epicentre is determined by the angular mo- 
mentum L of the star and the radial amplitude by the energy 
excess AE. 

The amplitude in tangential direction is 

2 



J/m = RoAlfn 



Po 



(8) 



The epicyclic ellipse is elongated into the radial direction, 
because 1^/3^2 for reasonable rotation curves. 



2.2 Mass loss of the star cluster 

We investigate the mass loss of a star cluster on a circu- 
lar orbit with Rc,Qc in the tidal field of the Galaxy. The 
orbit of the cluster is offset to the epicentres Rq , Qq of the 
stars in the tidal tails. The differential rotation results in 
the elongation of the tidal tails. 

As long as the stars are influenced by the gravitational 
potential of the cluster, it is appropriate to use a reference 
frame corotating with the star cluster. We use cylindrical 
coordinates (R, <p) corotating with the star cluster and with 
the origin at the galactic centre. The angular speed is Qc and 
the star cluster centre is at (Rc,^ = 0). Then we shift the 
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Figure 2. Top: Effective potential of a star cluster in the corotat- 
ing frame. LI and L2 are the Lagrange points. Bottom: Sketch of 
escaping stars in the effective potential of a star cluster through 
the Lagrange points LI and L2. 



origin to the star cluster centre and use a local cylindrical 
coordinate system with x — R — Rc and y = Rip (for an 
illustration see figure [T]). 

In the corotating system we set the zero point of 
the galactic potential $ 9 such that the effective potential 
* E ,eff(JJc) = $ B - Q. 2 c Rq/2 = vanishes at R c . We get 



o 2 

"C / n 2 r,2 



(9) 



The Jac obi energy Ej, which is th e only known constant of 
motion ijBinnev fc Tremaindll987h . is given by 



Ej = E - n c L = $ cff + 



(10) 



with velocity v in the corotating rest frame. The Jacobi en- 
ergy of the cluster motion in the galactic potential vanishs. 
The effective potential has saddle points at the inner and 
outer Lagrange points L1/L2 (Fig. [2}, where the stars are 
leaking out. All stars with Jacobi energy exceeding the crit- 
ical value -Ej,crit = & c ft(Ll/L2) can in principle leave the 
cluster. 

For the effective potential we get at y = 



$ c fr(x,0) 



2 



ficz +*d(a:,0) 



(11) 



The tidal radius j"l is given by the distance of the Lagrange 



points to the cluster centre. It is determined by 



= 



d$„ 



da; 



"c(/%-4)r L ± 



GA/ cl (r L ) 



where M c \(rt) is the cluster mass enclosed in 
the well known equation for the tidal radius 

±GM cl 



(12) 



We find 



(13) 



(4 ~ H C J 

where we assumed that the full cluster mass is enclosed in 
The effective potential at the Lagrange points is 



$ e ff(|r L |,0) = Ej,. 



-0c)^c r t = - 



3 GM C 



(14) 



2 |fx| 

The last expression shows that the contribution from the 
star cluster potential is twice that of the effective potential 
of the Galaxy. 

A star starting near LI or L2 with velocity dl escapes at 
constant Jacobi energy but with changing energy and angu- 
lar momentum until the cluster potential can be neglected. 
Then the position (x, y) and velocity v — (v r ,Vt) in the tidal 
tail are related to (rL,0) and vl by 



£j = $ cff (r L ,0) 



leading to 



3^ + 



= 3 + 



A(v 2 



f(*) 



V 

Y 



(4-/3g)fi 2 c 



GM c i/|r L | 



(15) 



(16) 



(17) 



■«L. 



Stars moving along the equipotential 
surface (A(v 2 ) — 0) yield as initial position x = \/3Vl and 
initial velocity essentially tangential vt ~ «l- This approxi- 
mation fits well with the radial position of the equipotential 
surface through L1/L2 at large distances from the cluster in 
Fig. [21 Stars moving radially gain kinetic energy (Aw 2 > 0) 
resulting in a larger x and stars starting tangentially loose 
kinetic energy (Av 2 < 0) leading to a smaller x. 

For a continuous mass loss until dissolution it is neces- 
sary that the Jacobi energy of bound stars is lifted above the 
critical value £j lCr it, which increases due to the mass loss. 
There are two physical effects, which are responsible for a 
continuous mass loss of the cluster. The first one is triggered 
by the mass loss of the cluster itself. Mass loss on a timescale 
large compared to the dynamical time of the cluster leads 
to an increase of Ej of the bound stars by 

But the critical value Ej lC tlt increases more slowly, because 
the tidal radius decreases with decreasing mass 

d$efr(rL,0) 



df 



oc M 



2/3 



(19) 



Initiated by mass loss due to stellar evolution or by a few 
stars above Ej lCI it mass loss will continue by stars lifted 
above the critical value. 

The second process is dynamical evolution of the clus- 
ter due to 2-body encounters. With the relaxation timescale 
stars are scattered above £Jj, cr it and can leave the cluster. 
The relative importance of the two effects depend on the 
mass, number of stars and the structure of the cluster. 
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2.3 Dynamic parameters of tidal tail stars 

Since the orbits are epicycles perturbed by the acceleration 
of the cluster, the connection of the initial position and ve- 
locity (rL, Vl) to (x,v) at a later time, when the cluster 
potential can be neglected, is very complicated. Here we are 
interested in the statistics of initial and final properties of 
the escaping stars. 

For the transition from bound stars to escaped stars, 
we need to combine the motion in the frame corotating 
with the cluster Rc,Qc and that in the non-rotating refer- 
ence frame, where we derived the properties of the epicycles 
around Ro,^lo- For measuring the shape and kinematics of 
the tidal tails we stay in the corotating rest frame centered 
at the cluster. Therefore we transform the epicyclic motion 
to the corotating frame with respect to Rc, fic- 

The radial offset ARo — Ro — Rc of the epicentre of 
a star is determined by the angular momentum difference 
AL = L — Lc (see equation I A9[) . Here we need only the first 
order term of ARo in AL, which is 



ARp 
Rc 



2 AL 



2x 



— ) 

Me) 



(20) 



Since the epicycles are counterrotating with respect to the 
disc rotation, the relative velocity in the tidal tails is smallest 
at the pericentres (with respect to the cluster motion) . These 
are the locations where the clumps occur. The tangential 
distance yo(T) of the pericentres in the corotating frame are 
determined by the shear flow of the epicentre motion. The 
period is T = 2-k/kq leading to 



2tt 



Ro (fio - fic) 



re 
2tt 0c 



2k Rc 

Pc 



13c S 
47vf3 2 c 



Pc 



31 



x + 



Pc 

Vt 

2Q C 



PI 



Lc 



(21) 
(22) 



For the special c ase of x = a nd y± — we recover the 
equation given in iKiipper et al\ ^2008). The spread in AL 
of the escaping stars leads to a corresponding tangential 
width of the first clump. In the succeeding clumps the tan- 
gential spread increases linearly and quickly smears out the 
pericentre positions over a whole period of yo(T). 

The radial spread of the tidal tail stars is determined 
by the combined spread in ARo and in the amplitudes 
r m . The amplitudes are determined by the second order 
terms in the energy excess of the stars with respect to the 
epicentre energy Eo- Relative to the cluster center energy 
E c = $ g (Rc) + 0. c R c /2 and with the help of equation |5] 
we get for the energy difference of the epicentre 



AEo 



Eo -E c = n c AL - 



p c nc AL i 



2Pl L C 



(23) 



The amplitude of the epicycle is determined by the en- 
ergy excess with respect to the epicentre energy Eo and can 
be calculated using the Jacobi energy of the star 



E.i = •f>g,e«(Rc+x) + Y 



- 0Cr<2 



Qfof + — (24) 



(taking into account the zero-point $ g , e ff(.Rc) = 0) and 



equation 1231 leading to 



AE 



E- E = Ej + n c AL - AEo 



Ej + Q C L C 



4_Pc AL 2 
2/3, 



T 2 



= E 3 + ^^^ARl 
Now the epicyclic amplitude r m is determined by 



2AE 

pm 



Pc a c>2 , 2 Ej 



(25) 
(26) 

(27) 
(28) 



Pi 



AL 2 | Ej 

L 2 C Q c 



The position of the apo/pericentre with respect to the clus- 
ter centre is given by 



Rn 



Rc = A7?o ± r n 



(29) 



In terms of position x and velocity v T ,Vt in the corotating 
frame the apo- and pericentre are given by 



1 / 2v t 

-m ( 4x + FT 

Pc V 



(30) 



± /T /((4 - 



Pl) 



2v_t 

Q.C/ " C 

If the star starts at peri/apo-centre (with no radial velocity 
v T = 0) we find the corresponding apo/peri-centre by 



^ r = : 



4(4* + &±|(4-^)*-&|) 



(31) 



In lKiipper et al\ (|2008l ) the special case of v r = vt = and 
x = tl was adopted, where the star has a Jacobi energy of 
E 3 = -GM c i/(2r h ) = Ej.crit/3. 



3 NUMERICAL MODELING 

For the main investigation we used the direct 0GRAPE N- 
body code to calculate the evolution of star clusters in an 
analytic Galaxy model. It is described in this section. For the 
case of a star cluster near the Galactic centre and for testing 
the reliability of the simulations we used a variant of the 
direct iV-body6++ code, which is described in Section \4. 41 



3.1 0GRAPE TV-body code 

For the high resolution direct TV-body simulations at large 
galactocentric distances we used the specially developed 
^GRAPE code. The code itself and also the special 
GRAP E hardware is described in more detail in lHarfst et al\ 
(2007). Here we mention briefly the most important spe- 
cial features of the code. The program was already well 
tested with different A^-body applications including the high 
resolution study of the dynamical evolution of the galactic 
centre with a binary (or s i ngle) Super-massive Black Hole 
jBerczik et al\ \200B . 120061 ; iMerritt et ai] |2007| ). The same 
code was also recently used to study the s hape parameters of 
a larg e set of rotating open star clusters dKharchenko et ali 
2008B 



1 The present version of the code will be pub- 
licly available from one of the authors FTP site: 
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Table 1. The list of galaxy component parameters. The first 
column gives the component, the second the mass, and the third 
and fourth the Plummer-Kuzmin parameters (equation [32} . 



Mass component 


M [M ] 


a [kpc] 


6 [kpc] 


Bulge 


1.4 ■ 10 10 


0.0 


0.3 


Disk 


9.0 ■ 10 10 


3.3 


0.3 


Halo 


7.0 • 10 11 


0.0 


25.0 



The program acronym <?!>GRAPE means: Parallel 
Hermite Integration with GRAPE. The serial and paral- 
lel version of the program has been written from scratch in 
ANSI-C and uses the standard MPI library for communica- 
tion. For the calculation of the star cluster dynamics in the 
galactic potential we use the parallel GRAPE systems built 
at the Astronomisches Rechen-Institut in Heidelberg, and 
at the Main Astronomical Observatory in Kie\{f|. 

The program uses the 4-th order Hermite integration 
scheme for the particles with hierarchical individual block 
timesteps, together with the parallel usage of GRAPE6a 
cards for the hardware calculation of the acceleration a and 
the first time derivative of the acceleration a (this term is 
usually called 'jerk' in the A r -body community). 

For the simulation of star clusters in the tidal field 
of the Galaxy an analytic external potential is added. We 
use an axi-symmetric three component model, where bulge, 
disc and halo are describ ed by Plummer-Kuzmin models 
l|Mivamoto fc Nagai|[l975n with 



9{R,z) = - 



G-M 



(32) 



R 2 + (a + Vb 2 + z 2 



where b is a measure of the core radius and a a measure of 
the flattening. 

For the parameter s we use similar values as in 
iDouphole fc Colinl (| 19951) with slightly corrected masses (see 
Table [TJ to reproduce the observed Milky Way rotation 
curve in the solar neighbourhood. The rotation curve, the 
epicyclic frequency (3 — /t/fi (normalized to the orbital fre- 
quency Q) and the logarithmic derivative /?' = d/3/d In R are 
shown in Fig. [3] 



3.2 Initial conditions for the star cluster 

The s tar clusters ar e modeled star by star using a Salpeter 
IMF (|Salpeterlll9"55T ) in the mass range of 0.08 < m/M Q sC 
8. We include a simple mode l for s tellar evolution from 
Ivan den Hoek fc Groenewegenl (1 19971 ) and distribute the 
stellar mass loss uniformly over th e metallicity dependent 
stellar lifetimes l|Raiteri et aZ.|[l996l ). 

For the generation of the initial particle distribution and 
velocities we use a nonrotating King model with Wo = 6.0. 
In the A r -body code the physical quantities are normalized 



f tp: //f tp. ari .uni-heidelberg. de/staf f /berczik/ 
phi-GRAPE-cluster /code-paper/. 

2 GRACE: http://www.ari.uni-heidelberg.de/grace 

3 GOLOWOOD: http://www.mao.kiev.ua/golowood/eng 
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Figure 3. Top: Rotation curve of the Galaxy model and the 
contributions from bulge, disk and halo. Bottom: Epicyclic fre- 
quency parameter /3 and the logarithmic derivative j3' . The ver- 
tical dashed lines mark the orbits of the star clusters selected for 
the investigation. 



by G = Mci = 1 and total e nergy E TOT = -0.25 as first 
introduced by I Aarseth et aZj(|l974h . 

After generating the dimensionless parameters for the 
cluster data we set the different physical mass & half-mass 
radius parameters for the clusters. In general these two pa- 
rameters can be set independently, but in order to reduce 
the number of free parameters we decided to use a phys- 
ically motivated relation between initial mass and initial 
half-mass radius of the star clusters. For this purpose we 
use the extension of the well known mass vs. radius rela- 
tion observed in molecular clouds and clumps in our Milky 
Way. For relevant re ferences see t h e list of observationa l 
and theo r etical p apers lLarsod(ll98lT):ISolomon et al\ (Il987l); 
iMalonevI l|l99fj|) : iTheis fc Henslerl l|l993h : llnoue fc Kamaval 
(2000). Using such an approximation we can write down a 
scaling relation for our initial cluster size and mass. For fix- 
ing the size of the cluster we use the radius containing 60% 
of the cluster mass (because this radius is approximately in- 
dependent of rotation for future extension to rotating King 
models). We set 



r k « 100 • 



10 6 Mf7 







[PC] 



(33) 



For the used three physical masses M c i — 10 s , 5-10 s 
and 10 4 [Mq] we get the corresponding radii rk = 3.0, 7.0 
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and 10 [pc]. We do not adapt the 'tidal radius' ol the King 
models, where the density vanishes, to the tidal radius tl 
of the galactic field. The starting point of the star cluster 
orbit is determined by the position in the galactic disk (0.0, 
Rc, 0.0) with corresponding velocity (-Vc, 0.0, 0.0) added 
to each star. The star cluster is nonrotating and has exactly 
the angular momentum VcRc- If Vc is exactly the circular 
velocity at Rc, then the angular momentu m corresponds to 
Lg of the circular orbit at 7?c • In contrast iFukushige et al\ 
(2000) started with clusters in the corotating frame lead- 
ing to an additional spin of the star cluster. For practical 
reasons the initial velocity was slightly smaller then the cor- 
reponding circular speed at Rc by neglecting the decimal 
places. Therefore the star clusters in models 01-09 started 
at apo-centre and moved on an epicycle relative to the clus- 
ter epicentre motion. It turned out that the galactocentric 
distance variations of the cluster orbits are comparable to 
the tidal radii of the star clusters. Since the Jacobi energy is 
conserved only in a rest frame with constant angular speed, 
the origin of the coordinate system must be determined by 
the epicentre motion of the cluster and therefore the cluster 
centre is moving on an epicycle in that coordinate system. 
The main effect on the evolution of the tidal tails is an ad- 
ditional periodic force with the epicyclic frequency. In order 
to test the effect of this 'resonant forcing' we set up model 
10 which is on an exact circular orbit. The only difference to 
model 08 is the larger initial velocity by 0.3 km/s. The dif- 
ferences between these two models in the mass loss rate and 
in the position and strength of the tidal clumps are negligi- 
ble. Therefore we use model 10 as the fiducial model for the 
detailed investigations. The other models are used to inves- 
tigate the parameter dependences of the tidal tail structure. 
The cluster parameters of the models are listed in Table d] 



4 RESULTS 

We discuss in detail the properties of model 10 on an ex- 
act circular orbit. In this case the origin of the corotating 
coordinate system is at the cluster centre. Firstly we an- 
alyze the mass loss and the Jacobi energy distribution of 
the bound stars. Then we determine Jacobi energy, angular 
momentum and energy excess of the stars in the tidal tails 
and compare the predictions with the clump positions and 
widths. Then we discuss the parameter dependence of the 
tidal tail structure. Finally we present a numerical simula- 
tion of a star cluster near the galactic centre to demonstrate 
the generality of the theory. 

4.1 Cluster mass loss 

There is no unique definition of bound stars for star clus- 
ters in tidal fields. The main reason is that many stars with 
Jacobi energy exceeding the critical value remain for a long 
time in the vicinity of the cluster. For determining the mass 
loss rate and for visualisation we use a rather conservative 
measurement by using an energy criterion in the comoving, 
but not corotating, reference systerrQ. We assume, that all 

4 The video snapshots from all the simulations 
will be publicly available from the FTP site: 
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Figure 4. Jacobi energy distribution of dynamical cluster mem- 
bers and of stars inside a sphere with radius |tl | at t=471 Myr of 
model 10. Note that lower values of Ej are to the right. 

particles which have a negative relative energy in the cluster 
potential are still bound to the cluster: 

|£i GRA |d > E? m . 

An inspection of the particle distribution shows that this cri- 
terion coincides approximately with stars inside \vl\. On the 
basis of this criterium we create the list of particles denoted 
by 'dynamical' cluster members. In Fig. U the difference in 
the Jacobi energy distribution of dynamical cluster mem- 
bers and of stars inside a sphere with radius |rx| arc shown. 
There are only small deviations at the high energy end. 

In Figure [5] we present the mass evolution of the star 
clusters. In the upper panel the models at Rc = 8.5 kpc are 
shown. Mass loss of model 08 at a slightly eccentric orbit 
shows a very small modulation compared to the correspond- 
ing model 10 on the exact circular orbit. The top line shows 
the mass loss of all stars due to stellar evolution for model 
10. The bottom panel shows the mass evolution of models 
with initial mass M c i = 10 4 Mq at different galactocentric 
distances. 

Figure [5] shows the effective potential of bound stars in 
units of initial GM a i/r^ at different times as a function of ra- 
dial position x covering a mass loss range of 50% to 90%. At 
the last timestep the cluster mass is already reduced to 10% 
of the initial mass. The full lines show the approximation of 
$eff(£,0, 0) from equation 1111 using a point mass potential 
for the star cluster. It is a lower boundary of $ e fr(a:, y, z) in 
projection and shows a perfect agreement in the vicinity of 
the tidal radius tl, which is marked by the tics. According 
to equations [13] and [14] the tidal radius rx and the critical 
energy -Ej, C rit scale with M^/ 3 (t) and M^/ 3 (t), respectively. 

Figure [7] shows the Jacobi energy Ej in units of initial 
GM c i/rL of all bound stars as in Figure [5] but as function of 
distance to the cluster centre rdc- The full lines are the same 
functions as in Figure [7] At all times there is a considerable 
number of stars exceeding the critical value -Ej, C rit which is 
marked by the dashed horizontal lines. These potential es- 
capers are well distributed all over the cluster. In Figure [5] 

ftp : //ftp . ari . uni-heidelberg . de/staf f /berczik/ 
phi-GRAPE-cluster/video-paper/pos/. 
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Table 2. The model parameters of all runs. Column 1 gives the number of the model, columns 2-4 are initial mass, number of particles 
and scale radius of the cluster, columns 5 and 6 are the initial position and velocity of the cluster, column 7 is the epicycle parameter at 
distance Rq, column 8 is the initial Lagrange radius of the cluster. 



# 


M o] [M ] 


N 


ffc [pc] 


R C [kpc] 


V C [km/s] 


P 


rL [pc] 


01 


10 3 


4040 


3.0 


7.0 


236 


1.363 


12.08 


02 


10 3 


4040 


3.0 


8.5 


233 


1.372 


13.93 


03 


10 3 


4040 


3.0 


10.0 


231 


1.396 


15.78 


04 


5 10 3 


20202 


7.0 


7.0 


236 


1.363 


20.66 


05 


5-10 3 


20202 


7.0 


8.5 


233 


1.372 


23.82 


06 


5-10 3 


20202 


7.0 


10.0 


231 


1.396 


26.98 


07 


10 4 


40404 


10.0 


7.0 


236 


1.363 


26.04 


08 


10 4 


40404 


10.0 


8.5 


233 


1.372 


30.00 


09 


10 4 


40404 


10.0 


10.0 


231 


1.396 


34.00 


10 


10 4 


40404 


10.0 


8.5 


233.297 


1.372 


30.00 
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Mass loss due to stellar evolution 
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Figure 5. The time evolution of bound mass of the different mod- 
els. Top: All models at distance Rc = 8.5 kpc. The top-most line 
shows the mass loss due to stellar evolution of all stars. Bottom: 
Models with initial mass M c \ = 10 4 Mq at different distances to 
the galactic centre. 
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Figure 6. The effective potential of all bound particles for dif- 
ferent times of model 10. The full (red) lines show the analytic 
approximation according to equation 1111 The vertical tics show 
the position and energy of the Lagrange points tl . 



lower panel shows the cumulative distributions of Ej start- 
ing at the high energy end. The critical values are marked 
by crosses showing that at all times up to the dissolution the 
fraction of 'potential escapers' is 35%. In order to measure 
the importance of 2-body encounters on the evolution of Ej 
more detailed star-by-star investigations are necessary. This 
is beyond the scope of this paper. 



we quantify the energy distribution of the bound stars. The 
histograms of Ej in the upper panel of Figure [8] show the 
temporal evolution of bound stars, which demonstrates that 
the maximum of the distribution is near the critical energy 
with a large fraction of stars above the critical energy. The 



4.2 Tidal tail clumps 

The structure of the tidal tails are determined by the an- 
gular momentum offset AL and energy excess AE of the 
stars. All stars move on epicycles with angular frequency 
k. The epicentres are determined by AL (see equation I20|l 
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6130.1 [Myr] ; 
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Figure 7. The Ej distribution of all bound stars as function of 
distance to the cluster centre 7\j c defined by the density centre 
for different times of model 10. The full (red) lines are the same 
as in Figure [6] The dashed (cyan) lines mark the critical Jacobi 
energy. 



and the amplitudes by AE (see equation I29[). Clumps form 
at the pericentres, where the streaming velocity is minimal. 
Figure [9] shows the clumps in the tidal tails for two differ- 
ent times. The local density at the positi on of each star is 
determ ined by the neighbour criterion of ICasertano fc Hutl 
l| 19851 ) using 10 neighbours. It is colour coded in logarithmic 
scale. Note that the clump density does not decrease with 
decreasing cluster mass. The circle marks the tidal radius. 
The structure is similar in the trailing and leading arm as 
expected from the symmetry in the epicyclic approximation. 

The radial offset ARo of the epicentres of the tidal tail 
stars is predicted by equation [20] from the measured AL. In 
Figure \W\ the histograms of calculated ARo scaled to the 
corresponding tidal radii ri. are shown for different times. 
For the time resolution we have selected stars at tangential 
distances 200 pc< \y\ <400pc. We find that ARo is propor- 

1/3 

tional to n, which scales with MJ . In the next section we 
quantify the scaling. The corresponding distributions of Ja- 
cobi energies Ej are shown in the lower panel of Figure 1101 
Here Ej is normalized to the actual GM c \(t) / Vl(4) • 

The tangential positions of the clumps are multiples of 
yo(T), which are connected to AL (equation I22[) . Figures [TT1 
and [T2] show a comparison of the density distribution along 
the tidal tails and predicted histograms of 1,2,3 x yo(T) 
for the early and late time. The position and width of the 
clumps agree well in both plots and the second and third 
clump show some overlap as predicted. 

A similar comparison of the apo- and pericentre distri- 
bution in the radial coordinate of the numerical simulation 
also matches the spread of the analytic prediction. 

In Figure Q3] the peri- and apocentre positions of the 
selected stars are shown as function of ARo- There is a 
strong correlation between epicentre offset ARo and ampli- 
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Figure 8. Top: The Jacobi energy distribution of bound particles 
at different times of model 10. Bottom: The cumulative distribu- 
tion of bound particles as function of decreasing Jacobi energy for 
different times. The critical Jacobi energy is marked by crosses. 



tude, that is of AL and AE. The double dotted (black) line 
shows the epicentre position and the dotted (blue) lines show 
the pericentre position with v% = and the corresponding 
apocentre (cycloids in the corotating reference frame). The 
dot-dashed (orange) lines are for Ej = which determines 
the maximum amplitude for most stars. Only very few stars 
fall ou tside this limit. The orbit adopted by iKiipper et al\ 
( 2008) with (x, v) = (rh, 0) is marked by the (blue) crosses. 
It is a typical orbit with a slightly smaller epicentre offset 
ARo compared to the mean value. Therefore they underes- 
timated the distance of the clumps slightly as was already 
obvious from their simple iV-body simulation. 



4.3 Parameter variation 

In models 01-09 of Table[2]we vary the cluster mass and the 
galactocentric distance. Model 10 at an exact circular orbit 
differs in the evolution of model 08 only in a small modula- 
tion of the mass loss. Therefore we discuss here models 01-09 
only. In the Section [4.41 a model near the galactic centre is 
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Figure 9. Positions and local densities of the stars of model 10 at an early and late time with about 50% and 10% bound mass, 
respectively. 



discussed to show that the theory holds also in this extreme 



4.4 The Galactic centre case 



Already in the temporal evolution of model 10 we found 
that the tidal tail structure scales with the tidal radius tl- 
For testing the dependence of ARo and yo(T) on the cluster 
mass M c i and on the distance Rc to the Galactic centre we 
measure for all nine models the value of yo(T) numerically. 
We add the trailing and leading arm density distributions 
and determine yo{T) by fitting the positions of clumps 1 and 
2 to the value of 1 x yo (T) and 2 x yo (T) . The result is plotted 
in Figure [T4l The upper panel shows the mass dependence. 
Best power law fits for each Rc give a power law index of 
1/3 to better than 1%. We test the scaling of ARo with tl 
by defining a scaling factor A 



ARo = Ar L 

and calculating A from yo(T) using equation 1211 



P yo(T) 



tt4- fP 



TL 



(34) 



(35) 



The result for the nine models is shown in the lower panel 
of Figure [T4l The values of A are independent of yo(T), that 
is of Mcl and R c . The best fit value is A = 2.77 ± 0.02. 
A possible dependence on f3 cannot be tested here, because 
the f3 values are very similar. A further discussion is given 
in Section l4~4l 



Near the Galactic centre the tidal forces are much stronger 
and the size of star clusters relative to the distance to 
the Galactic centre is much larger. Therefore it is an in- 
teresting case to test the epicyclic approximation in this 
extreme regime. We use the code NBODY6GC to simu- 
late the evolution of a star cluster in the tidal field of 
the Galactic centre. This code is based on the parallel 
A-bo dy code nbody6++ jAarsethl 1 19991. 120031; ISpurzeml 
Il999l ) and in detail described in lErnst et all ||2008h . The 
orbits of the stars in the st ar cluster are followed with 
a 4th-order Hermite scheme (jMakino fc Aarseth|[l992l ) m- 
cludi ng Kustaanheimo-Stiefel regula rization of close encoun- 
ters (jKustaanhcimo & Stiefel 19651 ) and Chain regulariza- 
tion l|Mikkola fc Aarsethl Il998l ). In addition, the orbit of 
the star cluster in the analytic background potential of 
the galactic ce ntre is followed using an 8th-ord er compo- 
sition scheme ( Y oshida 199(j; iMcLachlanl [l995) including 
the Chandrasekhar dynamical friction force w ith a variable 
Coulomb logarithm (|just fc Pefiarrubiall2005h . The dissipa- 
tive force is nu merically implemented wi th an implicit mid- 
point method l|Mikkola fc Aarseth|[2002T) . 

For the Galactic centre, we used a sca le free model 
e., with a power law density profile, e.g. iMezger et al. 



1996)) with a supermassive black hole l|Eisenhauer et al. 
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Table 3. The list of tidal clump parameters for the Galactic centre case as described in the text. 



Arm 


Clump 


AL/Lc 


<P0 [deg.] 


<Phi[deg.] 


¥>L2 [deg.] 


*<p/<P0 [%\ 


■r L [pc] 


yo [pc] 


Ay) 


-4l2 


AA/A L2 [%] 


Leading 


1 


-0.2179 


51.6 


43.3 


44.1 


14.5 


2.45 


16.9 


1.81 


1.54 


14.9 




2 


-0.2242 


55.3 


44.5 


45.4 


17.9 


2.76 


18.1 


1.72 


1.40 


18.6 


Trailing 


1 


0.2835 


65.8 


56.3 


54.8 


16.7 


2.45 


21.6 


2.31 


1.95 
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Figure 10. Top: The histograms show the distribution of the 
radial offsets N(ARq) °f model 10 at different evolution times 
normalized to the corresponding tidal radii rj^(t). Stars at dis- 
tances between 200 pc< \y\ <400pc along the tidal arms were 
selected. Bottom: The Jacobi energy distribution N(E]) of the 
same stars normalized to the actual GM c i(f)/ri,(t). Note that TV 
is the number of particles per bin and the bin width varies due 
to the different scaling. 



120051 ) added at the centre. The cumulative mass profile is 
given by M(R) = M hil + M (R/R ) a . The parameters were 
taken to roughly match those in the centre of the Milky 
Way. We used M hh = 3.6 x 10 6 M Q , M = 1.67 x 10 8 M Q 
at Ro — 20 pc and a — 1.2. At a distance of i?c=20pc 
the influence of the central black hole can be neglected and 
the model is scale free. In the limit of a scale free model, 
the ratio /3 = is independent of galactocentric distance 
and is given by /3| = a + 1 = 2 .2. For the star cluster, 
we used a King model ((King 1966) with Wo ~ 6, mass of 
M c i = 10 6 M Q , and half-mass radius of — 1.64 pc start- 
ing at a galactocentric radius of Rc = 20 pc. The particle 



Figure 11. Density and yo{T) plot along the tidal tails at 
t=1.4Gyr of model 10. Grey (green) dots are the local densi- 
ties and black (blue) dots show the mean density averaged over a 
30 pc sphere. The horizontal line corresponds to the stellar density 
in the solar neighbourhood. 



number is N = 10 5 . The initial tidal radius is tl = 2.72 pc 
(mean of LI and L2). 

Figure [TBI shows the tidal arms for the simulation after 
an evolution time of £=1.3 Myr. At £=1.3 Myr the galacto- 
centric distance has slowly decayed due to dynamical friction 
to Rc = 18.8 pc and the tidal radius decayed mainly due to 
cluster mass loss to tl = 2.26 pc according to equation [13] 
(both marked by circles). The local density is colour coded 
showing clearly the density maxima in the tidal tails. Three 
clumps can be identified in the leading arm and two clumps 
in the trailing arm. The radial lines from the Galactic centre 
mark the angles of three of these clumps with respect to the 
cluster centre. 

Figure [16] shows the density along the tidal arms as 
a function of azimuth angle tp with respect to the cluster 
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Figure 12. Same as in Figure ITT1 at time t=6.1 Gyr. 



centre. The clumps can be clearly identified as peaks. The 
^-range exceeds 360° showing the wrap of the tidal tails. 

Figure [T7] shows the histogram of the number of stars 
as function of angular momentum difference AL/Lc- We in- 
cluded only stars within an angle of 25 degrees around the 
density maxima in the statistics. The distribution is asym- 
metric with respect to the leading and trailing arms which 
cannot be explained in the frame of the epicyclic theory. The 
maxima in the histogram correspond to the density maxima 
of the clumps in the leading and trailing arms, respectively. 

The measured angles ipo = yo{T)/Rc of the density 
maxima in the clumps can be compared to the theoretical 
estimates ipL detrived from AL. Note that these are the 
angles between the cluster centre and the first clump or the 
first and the second clump. For the theoretical estimate, we 
plugged the most frequent angular momentum differences of 
the leading and trailing arms from Figure [17] into equation 



2n 

0s 



1 - 



Pi—4: AL 



Pi - 2 (3£ - 4 AL 2 



0i 



Lc 



L 2 c 



(36) 
(37) 



Compared to equation [2U we added here the second order 
term to test the sensitivity of the results. The results are 
shown in Table [3] for the first two clumps in the leading 
arm and the first clump in the trailing arm. The density 
maximum of the second clump in the trailing arm is not 
well-defined. We denoted the measured angle as cpo. The 
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Figure 13. Apo- and pericentre positions of tidal tail stars of 
model 10 as function of Ai?o. The double-dotted (black) line 
shows the epicentre position, the dotted (blue) lines are the peri- 
and apocentre positions with zero velocity at pericentre in the 
corotating frame, and the dot-dashed (orange) lines are peri- and 
apocentre positions wi th _Ej = 0. The cro ss marks the orbital 
parameters adopted in lKiipper et al ] l|2008h . 
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Figure 14. Top: Mass dependence of yo(T) for the nine different 
models. The lines show power law best fits of i/q. The points 
and lines are sampled by galactocentric distance. Bottom: Scaling 
factor A = ARqAl f° r models 01-09 with the same coding as 
above and the best fit value for A. 



theoretical estimates from the measured AL/Lc have been 
denoted as ip^i (first-order) and <£l2 (second-order). We de- 
fined Aip = ipo — ip-L2 for the error. Note that ipL2 is sys- 
tematically lower than ipo, but in principle we find a good 
agreement between measurement and theory. We also calcu- 
late the A values ^4l2 and A v from equations (|34[) and (1351) . 
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Figure 15. Clumps in the tidal tails of a star cluster with N = 
10° particles near the Galactic centre at t = 1.3 Myr. The color 
coding shows the local stellar density in logarithmic scale. The 
solid lines show the tidal radius and assigns the galactocentric 
distance of the cluster. The radial lines mark the angles of the 
first and second clump of the leading arm and the first clump of 
the trailing arm with respect to the cluster centre. 
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Figure 16. Density along the tidal arms as a function of the 
angle with respect to the cluster centre. The density maxima can 
be clearly identified. Stars inside the tidal radius corresponding 
to bound stars are in black (blue), the leading arm in dark grey at 
positive angles (red) , and the trailing arm in light grey at negative 
angles (green). The wrap of the tidal tails can be seen from the 
1,3-range exceeding 360°. 
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Figure 17. Histogram of angular momentum differences AL/Lq 
scaled by the angular momentum of the circular orbit. JV is the 
number of particles with a certain angular momentum difference. 
Only stars within 25 degrees around the density maxima in the 
clumps have been included in the statistics. 



^4l2 is calculated from the Taylor expansion of ARq in L to 
second-order. For A v we used the measured j/o = Rc<Po- We 
used the tidal radius according to equation (|13|) at the time 
t — cT, where the stars in the clumps were released from the 
cluster. Here T is the epicyclic period at Rc and c = 1 and 
c — 2 correspond to the first and second clump, respectively. 
We defined AA = — Al2 for the error. Again, the agree- 
ment is relatively good. However, the A value is significantly 
smaller than that at large galactocentric distances. 

There are some aspects of the tidal tail structure which 
cannot be explained by the simple epicyclic theory. The most 
prominent one is the asymmetry in the tidal arms concern- 
ing the angular momentum and Jacobi energy distribution. 
The errors Atp/tpo and AA/Al2 stem from a slight non- 
conservation of angular momentum in the tidal arms. The 
reason is the influence of the cluster potential. This needs 
to be investigated further. The asymmetry between the in- 
ner and outer Lagrange points with respect to the central 
potential of the cluster is only a few percent. However, due 
to the proximity to the Galactic centre, the phase space for 
the particles which escape into the leading arm is consid- 
erably smaller than for those which escape into the trailing 
arm. The streaming velocitiy differs considerably between 
the leading and trailing arms. Thus the redistribution of en- 
ergy and angular momentum for fixed Ej can be different. 
Since the radial offset is not small compared to the distance 
to the Galactic centre, an epicyclic theory for larger am- 
plitudes would be helpful. A Taylor expansion in rj = 1/R 
holds up to ecc entricities of 0.5 as was s hown by iDekkerl 
i|l97fj ) (see also lArifvanto fc Fuchsl (|2006h for a derivation 
in the s olar neighbourhood ). For a further investigation we 
refer to lErnst et al\ l|200ct ). where the effects due to third- 
and higher-order terms in the Taylor expansion of the effec- 
tive potential are discussed. 



5 SUMMARY 

We presented a quantitative derivation of the angular mo- 
mentum and energy distribution of escaping stars from a 



Clumps in the tidal tails of star clusters 13 



star cluster in the tidal field of the Milky Way. Despite the 
motion on a circular orbit, the tidal tails are clumpy due to 
the epicyclic motion of the stars. We compared the derived 
distances and widths of the clumps with numerical simula- 
tions using star-by-star simulations. For star clusters at the 
solar circle we included an IMF and mass loss due to stel- 
lar evolution in the calculations. The same equations were 
applied to a star cluster very close to the Galactic centre, 
where the tidal forces are very strong. 

We find a very good agreement of theory and mod- 
els concerning the tidal tail structure. The positions of the 
clumps are determined by the angular momentum offset 
of the stars, which lead to a radial offset of the epicen- 
ters with respect to the cluster orbit. The investigation of 
iKiipper et al\ |2008) is a special case of our investigations 
but for a Kepler potential. We find that the radial offset 
of the tidal arms is proportional to the tidal radius. How- 
ever near the Galactic centre the factor of proportionality is 
considerably smaller. 

The tidal arm structure at large galactocentric radii is 
symmetric, whereas the asymmetry near the galactic centre 
is considerable. This can be reproduced only partly by taking 
into account the correction of the epicyclic frequency at the 
epicentre radii. 

We have also measured the Jacobi energy distribution of 
bound stars and showed that there are 35% of stars above 
the critical Jacobi energy independent of the evolutionary 
state of the cluster. These stars can potentially leave the 
cluster. This is a hint, that mass loss is dominated by a 
self-regulating process of increasing Jacobi energy due to 
the diminishing gravitational potential of the star cluster 
induced by the mass loss itself. 

Finally we consider the observability of the predicted 
clump properties in the tidal tails of star clusters. The iden- 
tification of tidal tail stars of open clusters on a circular or- 
bit is strongly hampered by the large number of nearby field 
stars with similar properties. But with differential methods 
using high quality data for distances and velocities it may 
be possible to identify the most prominent first clumps. We 
have shown that the maximum density in the first clump 
does not decrease with time until dissolution of the cluster. 
The first clumps are formed by escaping stars with a time 
delay determined by the epicyclic period T K m 150 Myr. 
The tidal tail structure will probably survive the gravita- 
tional scattering proc ess, which is also responsible for the 
galactic disc heating (|Wielenlll977l ). On a timescale of one 
epicyclic period, we expect only small perturbations of the 
tidal clump position and velocity but no destruction. The 
density maximum of the first clump is of the order of a few 
percent of the field density of the galactic disc. The velocity 
imprint by the epicyclic motion is of the order of 2 km/s. An 
overdensity with these properties may be difficult to observe 
for clusters on exact circular orbits, if there is no additional 
separating property. For young star clusters the high frac- 
tion of early type stars can serve for such a discrimination. 
On the other hand most star clusters are identified by the 
systematic peculiar motion with respect to the field stars. If 
the positions and velocities of the tidal clump stars are prop- 
erly predicted, they may be observable as moving groups. 

For an application of the tidal clump theory to the 
eccentric orbits of globular clusters a perturbation theory 
with respect to the 'free falling' comoving coordinate sys- 



tem would be necessary. In this case the Jacobi energy is no 
longer a constant of motion. This will be a matter of future 
investigations. Some numerical test runs have shown that 
tidal tail clumps are formed also on highly eccentric orbits 
without an additional external perturbation, but the geom- 
etry and density is modulated along the orbital position. 
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APPENDIX A: TAYLOR EXPANSIONS 

We use Taylor expansions of the radial variation of different 
quantities in the gravitational field of the Galaxy in R and 
L with respect to some circular orbit with Ro and Lo- The 
potential $ c i of the star cluster is not expanded. The Taylor 
expansions are applied to: 1) the energy of stars on eccen- 
tric orbits with fixed L in the galactic field; 2) the effective 
potential in the cluster frame; 3) the epicentre position and 
energy excess of tidal tail stars in the cluster frame. For the 
Taylor expansions we use 
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with angular speed Q and epicyclic frequency k. 

For the epicyclic motion in the tidal tails fSec. 12. H and 
I4.4|l the gravitational potential of the galaxy is needed in 



terms of r = R — Ro- We find to third order 
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Ro + i" m we find for 



For the apo- and pericentre R 
the kinetic energy 
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leading to the energy excess AE — E — Eq (relative to the 
circular orbit with the same angular momentum Lo) 
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For the derivation in the rest frame of the star cluster 
we need the Taylor expansions of R(L), il(L) and ( l > g(L) 
with respect to AL = L — Lc (L is the angular momentum 
of the circular orbit at radius R(L)). In terms of position x 
and velocity v T ,v t in the corotating frame we have 

AL = (v t + fl c R)R~ «ci?c = {2tt c x + v t )R c + fici 2 + v t x 



The Taylor expansions in L axe 
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Here already second order terms contain f3' . Since 00' = 
d 2 In fi/d(ln -R) 2 , the logarithmic derivative (3' vanishes only, 
if Q(R) oc R a is exactly a power law with constant a. For 
realistic rotation curves varies considerably (see Fig. [3}. 
We find the energy excess AEq — Eq — Ec of the epicentre 
motion for equation 1231 by adding equations IA10I and IA11I 

1 AT 2 
AE = QcAL + — 2- (0c — 4) Clc—f — (A13) 

Here the (5' term vanishes. 

For the derivation in the rest frame of the star cluster we 
need $ e ff only in terms of r = R — Rc- The cluster potential 
is not expanded in a Taylor series. We find 

$es(R) = $el + * K ,«ff(Jfc) + — e 
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